
Bienvenides a la cuarta tarea del curso Statistical Thinking. Esta tarea tiene como objetivo evaluar los contenidos teóricos de la primera parte del curso, los cuales se enfocan principalmente en introducirlos en la estadística bayesiana. Si aún no han visto las clases, se recomienda visitar los enlaces de las referencias.
La tarea consta de una parte práctica con el fin de introducirlos a la programación en R enfocada en el análisis estadístico de datos.
Slides de las clases:
Videos de las clases:
Documentación:
En la siguiente sección deberá resolver cada uno de los experimentos computacionales a través de la programación en R. Para esto se le aconseja que cree funciones en R, ya que le facilitará la ejecución de gran parte de lo solicitado.
Para el desarrollo preste mucha atención en los enunciados, ya que se le solicitará la implementación de métodos sin uso de funciones predefinidas. Por otro lado, Las librerías permitidas para desarrollar de la tarea 4 son las siguientes:
# Manipulación de estructuras
library(tidyverse)
## ── Attaching core tidyverse packages ──────────────────────── tidyverse 2.0.0 ──
## ✔ dplyr 1.1.4 ✔ readr 2.1.5
## ✔ forcats 1.0.0 ✔ stringr 1.5.1
## ✔ ggplot2 3.5.1 ✔ tibble 3.2.1
## ✔ lubridate 1.9.3 ✔ tidyr 1.3.1
## ✔ purrr 1.0.2
## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## ✖ dplyr::filter() masks stats::filter()
## ✖ dplyr::lag() masks stats::lag()
## ℹ Use the conflicted package (<http://conflicted.r-lib.org/>) to force all conflicts to become errors
library(dplyr)
library(tidyr)
library(purrr)
# Para realizar plots
library(scatterplot3d)
library(ggplot2)
library(plotly)
##
## Adjuntando el paquete: 'plotly'
##
## The following object is masked from 'package:ggplot2':
##
## last_plot
##
## The following object is masked from 'package:stats':
##
## filter
##
## The following object is masked from 'package:graphics':
##
## layout
# Manipulación de varios plots en una imagen.
library(gridExtra)
##
## Adjuntando el paquete: 'gridExtra'
##
## The following object is masked from 'package:dplyr':
##
## combine
# Análisis bayesiano
library("rethinking")
## Cargando paquete requerido: rstan
## Warning: package 'rstan' was built under R version 4.4.2
## Cargando paquete requerido: StanHeaders
## Warning: package 'StanHeaders' was built under R version 4.4.2
##
## rstan version 2.32.6 (Stan version 2.32.2)
##
## For execution on a local, multicore CPU with excess RAM we recommend calling
## options(mc.cores = parallel::detectCores()).
## To avoid recompilation of unchanged Stan programs, we recommend calling
## rstan_options(auto_write = TRUE)
## For within-chain threading using `reduce_sum()` or `map_rect()` Stan functions,
## change `threads_per_chain` option:
## rstan_options(threads_per_chain = 1)
##
## Do not specify '-march=native' in 'LOCAL_CPPFLAGS' or a Makevars file
##
## Adjuntando el paquete: 'rstan'
##
## The following object is masked from 'package:tidyr':
##
## extract
##
## Cargando paquete requerido: parallel
## Cargando paquete requerido: dagitty
## rethinking (Version 2.01)
##
## Adjuntando el paquete: 'rethinking'
##
## The following object is masked from 'package:purrr':
##
## map
##
## The following object is masked from 'package:stats':
##
## rstudent
Si no tiene instalada la librería “rethinking”, ejecute las siguientes líneas de código para instalar la librería:
install.packages("rethinking")
## Warning: package 'rethinking' is in use and will not be installed
En caso de tener problemas al momento de instalar la librería con el código anterior, utilice las siguiente chunk:
install.packages("shape")
install.packages(c("mvtnorm","loo","coda"), repos="https://cloud.r-project.org/",dependencies=TRUE)
options(repos=c(getOption('repos'), rethinking='http://xcelab.net/R'))
install.packages('rethinking',type='source')
Un conjunto de carteros aburridos de las mordidas de perros ha
decidido realizar un catastro de mordidas recibidas por los empleados de
su empresa en un periodo de 8 meses, planeando en base a estos datos
realizar inferencia bayesiana. Los datos de las mordidas estas datos por
el dataset no+mordidas.csv, en donde cada fila representa
las mordidas recibidas por diferentes carteros y las columnas señalan si
fue mordido o no el cartero en los meses de estudio (si fue mordido es
señalado con un 1, de lo contrario es señalado con un 0). Cabe señalar
que un cartero no puede ser mordido mas de una vez al mes, ya que el
damnificado recibe licencia de todo un mes tras la mordida,
reincorporándose el siguiente mes al trabajo.
dataMordidas <- read.csv("no_mordidas.csv", header = TRUE)
cols <- c("bites_month_1", "bites_month_2", "bites_month_3",
"bites_month_4", "bites_month_5", "bites_month_6",
"bites_month_7", "bites_month_8")
grid_len <- 100
p_grid <- seq( from=0, to=1, length.out= grid_len)
prior <- dnorm(p_grid, 0.6, 0.05)
posteriors <- c()
counted.months <- c()
p_grids <- rep(p_grid, 4)
for (n in c(1, 2, 4, 8)) {
x <- sum(dataMordidas[, cols[1:n]])# total de mordidas
likelihood <- dbinom(x, n * nrow(dataMordidas), p_grid)
unstd.posterior <- likelihood *prior
posterior <- unstd.posterior / sum(unstd.posterior)
posteriors <- c(posteriors, posterior)
counted.months <- c(counted.months, rep(n, grid_len))
}
post.df <- data.frame("grid"=p_grids, "posterior"=posteriors, "months"=counted.months)
ggplot(post.df) +
geom_line(aes(x=grid, y=posterior)) +
facet_grid(months ~ .)
Respuesta:
Las curvas van adquiriendo una mayor altura y ocurriendo en un rango menor de valores para el eje x, esto supone que los datos que son usados para el modelo se van concentrando cada vez más en un rango más acotado de valores.
prior <- rep(1, grid_len)
posteriors <- c()
for (n in c(1, 2, 4, 8)) {
x <- sum(dataMordidas[, cols[1:n]])
likelihood <- dbinom(x, n * nrow(dataMordidas), p_grid)
unstd.posterior <- likelihood *prior
posterior <- unstd.posterior / sum(unstd.posterior)
posteriors <- c(posteriors, posterior)
}
post.df$posterior2 <- posteriors
ggplot(post.df) +
geom_line(aes(x=grid, y=posterior)) +
geom_line(aes(x=grid, y=posterior2), color="red") +
facet_grid(months ~ .)
Respuesta:
Las distribuciones para el prior uniforme están más a la izquierda que las de la prior gaussiana, esto ocurre ya que el prior uniforme asigna la misma probabilidad a todos lo datos, en cambio con el prior gaussiano se hace más probable el hecho de ser mordido por un perro y es por esto que se encuentra más a la derecha.
for (n in c(1, 2, 4, 8)) {
res <- quap(
alist(
W ~ dbinom (W+L, p) , # binomial likelihood
p ~ dunif(0,1)# uniform pior
),
data = list(W=sum(dataMordidas[, cols[1:n]]),L=n*250-sum(dataMordidas[, cols[1:n]])) )
df <- precis(res)
curve(dnorm(x, df$mean, df$sd), lty=2, add=FALSE)
}
Respuesta:
En este caso para todos los rangos de meses las distribuciones se ven centradas en un valor de 0.4 y al utilizar más datos para el modelo se van viendo menos casos lejos del centro, cada vez con “colas” más pequeñas, esto sucede ya que el metodo de laplace solo ve la punta que distribuye como normal y por tanto se tendrán valores más acotados, es decir la curva ocurrirá en rangos más “delgados”.
[] Grafique la densidad de la posterior y encuentre la proporción de los siguientes defined boundaries:
¿Cómo puede interpretar los resultados?
grid_len <- 1000
p_grid <- seq( from=0, to=1, length.out= grid_len)
prior <- rep(1, grid_len)
likelihood <- dbinom(6, 10 , p_grid)
unstd.posterior <- likelihood *prior
posterior <- unstd.posterior/sum(unstd.posterior)
samples <- sample(p_grid, prob=posterior, size=1e4, replace=TRUE)
dens(samples)
Intervalos:
# fracción de valores <0.35
sum(posterior [p_grid < 0.35])
## [1] 0.05002864
# fracción de valores entre 0.35 y 0.45
sum(posterior [0.35 < p_grid & p_grid < 0.45])
## [1] 0.1236835
# fracción de valores >0.45
sum(posterior [ 0.45 < p_grid])
## [1] 0.8262878
Respuesta:
Cerca de un 5% del área de probabilidad del posterior está bajo un
valor de 0.35.
Cerca de un 12,4% del área de probabilidad del posterior está entre un
valor de 0.35 y un valor de 0.45.
Cerca de un 82,6% del área de probabilidad del posterior está sobre un
valor de 0.45.
Lo cual sugiere que la concentración de los resultados está ligeramente inclinada hacia la derecha, tal y como se ve en el gráfico de la distribución de la posterior.
-[] Calcule un intervalo de credibilidad al \(50%\), \(75%\) y \(95%\). ¿Como se puede interpretar los resultados? ¿Cual podría ser un problema al usar intervalos de credibilidad?
Intervalo \(50\%\):
quantile(samples, c(0.25, 0.75))
## 25% 75%
## 0.4884885 0.6826827
PI(samples, 0.5)
## 25% 75%
## 0.4884885 0.6826827
Intervalo \(75\%\)
quantile(samples, c(0.125, 0.875))
## 12.5% 87.5%
## 0.4184184 0.7427427
PI(samples, 0.75)
## 12% 88%
## 0.4184184 0.7427427
Intervalo \(95\%\)
quantile(samples, c(0.025, 0.975))
## 2.5% 97.5%
## 0.3103103 0.8278529
PI(samples, 0.95)
## 3% 98%
## 0.3103103 0.8278529
Respuesta:
Hay un 50% de probabilidad de obtener un valor entre 0.49 y
0.68.
Hay un 75% de probabilidad de obtener un valor entre 0.42 y 0.74.
Hay un 95% de probabilidad de obtener un valor entre 0.3 y 0.83.
Es un rango de valores que creemos que el parámetro va a tomar con una cierta probabilidad.
Un problema puede surgir si las distribuciones obtenidas del posterior se comportan de una manera muy asimetrica, tal que las colas arrojen resultados engañosos y no congruentes a los intervalos requeridos.
Intervalo \(50\%\):
HPDI(samples, 0.5)
## |0.5 0.5|
## 0.4874875 0.6816817
Intervalo \(75\%\)
HPDI(samples, 0.75)
## |0.75 0.75|
## 0.4284284 0.7507508
Intervalo \(95\%\)
HPDI(samples, 0.95)
## |0.95 0.95|
## 0.3133133 0.8298298
Respuesta:
Dio como resultado valores bastante similares a los resultados de la parte anterior, esto era esperable ya que la distribución es bastante simétrica y por tanto las “colas” generadas en ambos serían bastante similares ya que la mayor concentración está en el centro.
los intervalos de HPDI buscan los intervalos más delgados y no necesariamento los que contengan la misma proporción de valores en las “colas”, además son más costosos.
# Número total de carteros y datos observados
num_carteros <- nrow(dataMordidas)
num_meses <- 8
total_mordidas <- sum(dataMordidas[, cols]) # total de mordidas
# Ratio real de mordidas
ratio.real <- total_mordidas
samples.sim <- rbinom(1e4, num_carteros*num_meses, prob= samples)
ggplot() +
geom_density(aes(samples.sim)) +
geom_vline(aes(xintercept=ratio.real), color="red")
Respuesta:
Este modelo no se ajusta bien a los datos, ya que la distribución de las predicciones no incluye la observación real como resultado central y probable, esto sucede ya que se tomó como hecho más probable el ser mordido, siendo la realidad que la frecuencia de carteros mordidos es mucho menor a lo asumido.
# Se escogen el mes 5 y 6
primer_mes <- dataMordidas$bites_month_5
segundo_mes <- dataMordidas$bites_month_6
# Carteros que fueron mordidos en el primer mes
mordidos_primer_mes <- which(primer_mes == 1)
cantidad_mordidos_primer_mes <- length(mordidos_primer_mes)
# Carteros mordidos en el segundo mes dado que fueron mordidos en el primero
mordidos_condicionados <- sum(segundo_mes[mordidos_primer_mes])
# Calculando la posterior con prior uniforme
grid_len <- 1000
p_grid <- seq( from=0, to=1, length.out= grid_len)
prior <- rep(1, grid_len)
likelihood <- dbinom(mordidos_condicionados, cantidad_mordidos_primer_mes , p_grid)
unstd.posterior <- likelihood *prior
posterior <- unstd.posterior/sum(unstd.posterior)
samples <- sample(p_grid, prob=posterior, size=1e4, replace=TRUE)
# Simulación de 10,000 réplicas
simulaciones <- rbinom(1e4, size = cantidad_mordidos_primer_mes, prob = samples)
ggplot() +
geom_density(aes(simulaciones)) +
geom_vline(aes(xintercept=mordidos_condicionados), color="red")
Respuesta:
Se puede ver una distribución bastante simétrica centrada en el valor del recuento real, por lo que se puede decir que son simulaciones bastante apegadas a la realidad de los datos del modelo.
En esta pregunta se trabajara con el dataset “notas.csv” el cual contiene las notas históricas de un curso desconocido. Suponga que los datos vienen de una distribución \(\mathcal{N}(\mu,\sigma^2)\), el objetivo de la pregunta es estudiar el comportamiento de los datos y los posibles valores de \(\mu,\sigma\) mediante técnicas de inferencia Bayesiana.
Usted sabe un dato extra sobre la información, \(\sigma \in [0.5,1.5]\) y, además, tiene una fuerte creencia a que es mas probable encontrar la desviación estándar real entre \([0.5,1]\) que en \((1,1.5]\).
# Leer información
data_notas <- read.csv("notas.csv")
head(data_notas)
# Función para crear likelihood dado mu y sigma
grid_function <- function(mu, sigma) {
# Log-likelihood basado en distribución normal
sum(dnorm(data_notas$Notas, mean = mu, sd = sigma, log = TRUE))
}
# Valores de la grilla
grid_mu <- seq(mean(data_notas$Notas, na.rm = TRUE) - 3,
mean(data_notas$Notas, na.rm = TRUE) + 3, length.out = 100)
grid_sigma <- seq(0.5, 1.5, length.out = 100)
# Crear grilla 2D
data_grid <- expand_grid(grid_mu, grid_sigma)
# Calcular likelihood para cada par (mu, sigma)
data_grid <- data_grid %>%
mutate(likelihood = map2_dbl(grid_mu, grid_sigma, grid_function))
# Priors
prior_mu <- rep(1 / length(grid_mu), length(grid_mu)) # Uniforme para mu
prior_sigma <- ifelse(grid_sigma <= 1, 0.75, 0.25) # Ponderación más alta para [0.5, 1]
prior_sigma <- prior_sigma / sum(prior_sigma) # Normalizar prior_sigma
# Crear columnas separadas para priors en data_grid
data_grid <- data_grid %>%
mutate(prior_mu = prior_mu[match(grid_mu, unique(grid_mu))],
prior_sigma = prior_sigma[match(grid_sigma, unique(grid_sigma))],
prior = prior_mu * prior_sigma)
# Calcular posterior no estandarizado
data_grid <- data_grid %>%
mutate(unstd_posterior = likelihood + log(prior), # Log-transformación para estabilidad numérica
posterior = exp(unstd_posterior - max(unstd_posterior))) %>% # Ajuste para evitar underflow
mutate(posterior = posterior / sum(posterior)) # Normalización del posterior
# Mostrar las primeras filas del dataframe resultante
head(data_grid)
Respuesta:
Prior para \(\mu\) (media):
Dado que no tenemos información previa sobre la ubicación de la media,
se asume una distribución uniforme para \(\mu\). Esto refleja una incertidumbre
inicial equitativa sobre los posibles valores de la media.
Prior para \(\sigma\) (desviación estándar):
Sabemos que \(\sigma \in [0.5, 1.5]\),
con una mayor probabilidad en \([0.5,
1]\) que en \((1, 1.5]\).
Para reflejar esta creencia, asignamos un prior con un peso de
0.75 para valores en \([0.5,
1]\) y 0.25 para valores en \((1, 1.5]\). Este prior sesgado refleja
nuestra expectativa sobre \(\sigma\).
# Punto de referencia
# Se recomienda cambiar estos valores por unos adecuados que le permitan estudiar
# Los valores de la distribución de mejor manera
valor_x <- 1
valor_y <- 1
# Grafico
punto_comparacion <- tibble(x = valor_x, y = valor_y)
plt <- data_grid %>%
ggplot(aes(x = grid_mu, y = grid_sigma)) +
geom_raster(aes(fill = posterior),
interpolate = T
)+
geom_point(x = valor_x, y = valor_y, size = 1.3,color="white")+
geom_label(
data = punto_comparacion, aes(x, y),
label = "Punto Comparación",
fill = "green",
color = "black",
nudge_y = 0, # Este parametro desplaza la caja por el eje y
nudge_x = 1 # Este parametro desplaza la caja por el eje x
)+
scale_fill_viridis_c() +
labs(
title = "Posterior para Mean y Standard Deviation",
x = expression(mu ["Mean"]),
y = expression(sigma ["Standar Deviation"])
) +
theme(panel.grid = element_blank())
plt
Respuesta:
Zona de máxima densidad: La mayor probabilidad para los parámetros \(\mu\) (media) y \(\sigma\) (desviación estándar) se concentra en un área definida. Esto indica que la estimación conjunta de los parámetros es consistente con valores específicos, sugiriendo un comportamiento estable en la media y la desviación estándar.
Forma de la distribución: La alta concentración de probabilidad en una región angosta sugiere baja incertidumbre sobre la estimación de \(\mu\) dentro del rango observado. La distribución está sesgada hacia un intervalo reducido para \(\sigma\), con la mayor densidad posterior centrada en valores cercanos a 0.75-1, lo que coincide con la creencia previa de que \(\sigma\) se encuentra en el rango \([0.5, 1]\).
# Codificamos los datos
x <- 1:length(data_grid$posterior)
# Sampleamos los indices
posterior_samples_aux <- sample(x,size = 1e4, replace = T, prob = data_grid$posterior)
# Obtenemos los verdaderos valores de la sampling distribution
posterior_samples <- data_grid[posterior_samples_aux,]
# Obtenemos solos los valores relevantes para la densidad
df <- data.frame(posterior_samples$grid_mu,posterior_samples$grid_sigma)
# Realizamos las densidades
dens(df)
Respuesta:
El proceso de sampling from grid approximation sirve para estudiar cómo se distribuyen los valores de los parámetros \(\mu\) y \(\sigma\) en la posterior. En el gráfico, se pueden observar las regiones de mayor densidad de puntos, lo cual sugiere que los valores de \(\mu\) y \(\sigma\) más consistentes con los datos se concentran en esas áreas. El gráfico muestra una mayor densidad en la parte baja del rango de \(\sigma\) (entre 0.55 y 0.7), esto confirmaría la prior que propusimos, que da mayor probabilidad a esa región. Además, la forma y dispersión de los puntos indican la certeza o incertidumbre en las estimaciones de los parámetros.
En cuanto a la distribución de la media (\(\mu\)), podemos observar que los valores más probables se encuentran en la región central del gráfico, donde se concentra la mayor densidad de puntos. Esto sugiere que la media estimada se encuentra dentro del rango 5.9 y 6.2, el cual indica la consistencia de los datos con los valores propuestos en la prior. Los puntos cercanos al centro de la grilla reflejan una mayor certeza sobre el valor de \(\mu\), mientras que aquellos más alejados de este centro reflejan mayor incertidumbre.
El objetivo de esta pregunta es introducirlo en la aplicación de una
regresión bayesiana. Con esto, buscaremos entender como calcular una
regresión bayesiana en base al “motor” aproximación de Laplace,
revisando como se comporta la credibilidad de sus predicciones y como la
regresión lineal puede llegar a mutar a aplicar una transformación en el
vector \(x\). Para responder esta
pregunta centre su desarrollo solo en las clases y las funciones que nos
brinda la librería rethinking.
Para esta pregunta usaremos el dataset Pokemon.csv, que
contiene las características de más de 700 personajes.
poke.df <- read.csv("Pokemon.csv")
head(poke.df)
Respuesta:
modelo <- alist(
Sp..Atk ~ dnorm(mu, sigma),
mu <- alpha + beta * Attack,
alpha ~ dnorm(100, 20), # Prior para alpha
beta ~ dnorm(0, 10), # Prior para beta
sigma ~ dexp(1) # Prior para sigma
)
El modelo propuesto busca modelar la relación entre
Sp..Atk (ataque especial) y Attack (ataque
físico). Las ecuaciones y los priors asumidos son los siguientes:
Modelo matemático:
\[ Sp..Atk_i \sim \mathcal{N}(\mu_i, \sigma) \quad \text{donde} \quad \mu_i = \alpha + \beta \cdot \text{Attack}_i \]
Priors definidos: - \(\alpha \sim \mathcal{N}(100, 20)\): Refleja
que esperamos un valor promedio de Sp..Atk alrededor de 100
con una incertidumbre razonable. - \(\beta
\sim \mathcal{N}(0, 10)\): Asume inicialmente que no hay una
relación fuerte entre Attack y Sp..Atk. -
\(\sigma \sim \text{Exponencial}(1)\):
Prior para la desviación estándar de los residuos, asegurando
positividad y favoreciendo valores pequeños.
fit <- map(
alist(
Sp..Atk ~ dnorm(mu, sigma),
mu <- alpha + beta * Attack,
alpha ~ dnorm(100, 20),
beta ~ dnorm(0, 10),
sigma ~ dexp(1)
),
data = poke.df
)
precis(fit)
El intercepto (\(\alpha\)) indica que, para un valor de Attack cercano a 0, el Sp..Atk predicho es 42.41. La pendiente (\(\beta\)) es 0.39, lo que implica que por cada aumento de una unidad en Attack, el Sp..Atk promedio aumenta en 0.39 unidades. El alto valor de \(\sigma\) (29.50) refleja una gran variabilidad no explicada por el modelo.
plot_data <- data.frame(Attack = seq(min(poke.df$Attack), max(poke.df$Attack), length.out = 100))
plot_data$mu <- coef(fit)["alpha"] + coef(fit)["beta"] * plot_data$Attack
ggplot() +
geom_point(aes(x = poke.df$Attack, y = poke.df$Sp..Atk), color = "red") +
geom_line(aes(x = plot_data$Attack, y = plot_data$mu), color = "orange") +
geom_ribbon(aes(x = plot_data$Attack, ymin = plot_data$mu - 1.96 * coef(fit)["sigma"],
ymax = plot_data$mu + 1.96 * coef(fit)["sigma"]), alpha = 0.3, fill = "orange")
Relación positiva moderada:
Attack y
Sp..Atk, lo que confirma que al aumentar el valor de
Attack, también aumenta el valor promedio de
Sp..Atk.Alta dispersión en los datos:
Attack) que influyen significativamente
en Sp..Atk.Intervalos de credibilidad amplios:
Valores atípicos:
Limitaciones del modelo lineal:
Attack y Sp..Atk. Un modelo más complejo
podría mejorar la capacidad predictiva.El modelo muestra que existe una relación positiva, pero no captura toda la variabilidad ni explica los datos de manera óptima debido a la alta dispersión y otros factores omitidos.
A work by CC6104